rm(list = ls())
warnings()
library(fOptions)

###########################################################
#II. MINIMAX STATISTIC----THEOREM 2----INDEPENDENT PROGRAM#
###########################################################
MM_function<-function(sigma, P){
  #MM_function is to calculate the difference between the dollar errors of two given options.
  #MM-function is used in the Minimax_stat_function.
  Call_mkt_1=unlist(P[1])
  Call_mkt_2=unlist(P[2])
  maturity_dt=unlist(P[3])
  X1=unlist(P[4])
  X2=unlist(P[5])
  S=unlist(P[6])
  risk_free=unlist(P[7])
  dividend=unlist(P[8]) 
    Call_BS_1=GBSOption(TypeFlag = "c", S = S, X = X1, Time =maturity_dt, r = risk_free,  
    b =risk_free, sigma = sigma)
    Call_BS_2=GBSOption(TypeFlag = "c", S = S, X = X2, Time =maturity_dt, r = risk_free,  
    b =risk_free, sigma = sigma)
    DEDiff=abs(Call_mkt_1 - Call_BS_1@price)-abs(Call_mkt_2-Call_BS_2@price)
  return(DEDiff)
}

##############################################
Minimax_stat_function <- function(sigma,ls_P) {  
  #Minimax_stat_function is to find the minimax statistic of two given options
  sigma_imp1=GBSVolatility(price=Call_mkt_1, TypeFlag="c", S=S, X=X1, Time=maturity_dt, 
  r=risk_free, b=risk_free, tol=0.0001, maxiter=100)
  sigma_imp2=GBSVolatility(price=Call_mkt_2, TypeFlag="c", S=S, X=X2, Time=maturity_dt, 
  r=risk_free, b=risk_free, tol=0.0001, maxiter=100)

  #Solve for sigma the function Dollar Error X1= Dollar Error X2, 
  #or MM_function =0 to find minimax statistic

print(sigma_imp1)
print(sigma_imp2)

  sigma_mm=uniroot(MM_function,P=ls_P, lower=min(sigma_imp1,sigma_imp2), upper=max(sigma_imp1,sigma_imp2),
  tol=0.000001)$root
  sigma=sigma_mm

  Call_BS_1_mm=GBSOption(TypeFlag = "c", S = S, X = X1, Time =maturity_dt, r = risk_free,  b =risk_free, sigma = sigma)
  #Call_BS_1_mm is the Black-Scholes call price of strike price X1 at sigma_mm

  Minimax_stat= abs(Call_mkt_1 - Call_BS_1_mm@price)
  return(Minimax_stat)
}

###############################
#REPORTING MINIMAX STATISTIC###
###############################
#Data reading and calculating Minimax statistics
d<-read.table('C:/Users/diepnguyen/Documents/MSF- Clark/Implied Binomial Tree/Research/R/OptionPairsMultiply.csv',header=TRUE,sep=",",na.string=".",as.is=TRUE)

v_minimax_stat<-vector(length = length(d$Call_mkt_1))
v_minimax_percentage<-vector(length = length(d$Call_mkt_1)) 
  #Minimax Percentage is the Minimax statistic divided by the stock price

for (i in 1:length(v_minimax_stat)) { #Run Minimax statistic function for each pair of options
  Call_mkt_1=d$Call_mkt_1[i]
  Call_mkt_2=d$Call_mkt_2[i]
  maturity_dt=d$maturity_dt[i]
  X1=d$X1[i]
  X2=d$X2[i]
  S=d$S[i]
  risk_free=d$risk_free[i]
  dividend=d$dividend[i]
  pair_category=str(d$pair_category[i])
  ls_ParamOptions=list(Call_mkt_1,Call_mkt_2, maturity_dt, X1,X2,S,risk_free,dividend) 

print(maturity_dt)
if (maturity_dt>0.1) {
  v_minimax_stat[i]=Minimax_stat_function(sigma,ls_ParamOptions)
  v_minimax_percentage[i]=100*v_minimax_stat[i]/S
}
else {
  v_minimax_stat[i]=0
  v_minimax_percentage[i]=100*v_minimax_stat[i]/S

}
}

#Printing the output to table and file
result_matrix <-data.frame(cbind(X1=d$X1, X2=d$X2, Call_mkt_1=d$Call_mkt_1, Call_mkt_2=d$Call_mkt_2,
  maturity_dt=d$maturity_dt, S=d$S, risk_free=d$risk_free,dividend=d$dividend,current_date=d$current_date,
  pair_category=d$pair_category, Minimax= v_minimax_stat, Minimax_percentage=v_minimax_percentage),
  stringsAsFactors = default.stringsAsFactors())

write.table(result_matrix, file = "C:/Users/diepnguyen/Documents/MSF- Clark/Implied Binomial Tree/Research/R/Result_Minimax.csv", 
  sep = ",", col.names = NA, qmethod = "double")

getwd() #See directory to the output file
